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Abstract. MHD turbulence has long been proposed as a mechanism for the heating of coronal loops in the framework 
of the Parker scenario for coronal heating. So far most of the studies have focused on its dynamical properties without 
considering its thermodynamical and radiative features, because of the very demanding computational requirements. In this 
paper we extend this previous research to the compressible regime, including an energy equation, by using HYPERION, a new 
parallelized, viscoresistive, three-dimensional compressible MHD code. HYPERION employs a Fourier collocation - finite 
difference spatial discretization, and uses a third-order Runge-Kutta temporal discretization. We show that the implementation 
of a thermal conduction parallel to the DC magnetic field induces a radiative emission concentrated at the boundaries, with 
properties similar to the chromosphere-transition region-corona system. 
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INTRODUCTION 

Magnetohydrodynamic turbulence in the framework of 
the Parker scenario for coronal heating lfTTl[T2l has been 
a very challenging problem to investigate numerically 
EH IS HUH [15]. As computers have advanced, it has 
become more feasible to do the large storage compress- 
ible problem. 

Why is it important to include compressibility and its 
related effects? There are three basic categories of inter- 
est with respect to active region loops and the coronal 
heating problem, viz.: structural, dynamical and thermo- 
dynamical. The most significant structural effect is strat- 
ification due to gravity. We can also modify this term to 
model the curvature of a typical loop. Among new dy- 
namical effects that are possible are compression and 
rarefaction of the plasma, as well as the formation of 
shocks. 

Thermodynamical effects include thermal conduction 
and radiation. In addition, the diffusivities can be tem- 
perature dependent. It's important to have these features 
in the model to begin to reproduce the energy cycle: ki- 
netic energy in the photosphere is transformed into mag- 
netic energy in the corona by means of photospheric foot- 
point convection. It is then transported to small scale by 
MHD turbulence, where through magnetic reconnection 
it is converted into thermal, kinetic and perturbed mag- 
netic energies. Heat is then conducted from the high tem- 
perature corona back toward the low temperature photo- 
sphere, where it is lost via optically thin radiation. 

Incompressible and cold plasma models only contain 
the first few parts of this energy cycle, without taking 
into account the thermodynamics. Any magnetic energy 



lost through Ohmic diffusion and any kinetic energy lost 
through viscous diffusion is simply lost from the system 
and the physics involved with thermal conduction and 
radiation is irrelevant. 

Our new compressible code HYPERION has allowed 
us to make a start at examining the fully compressible, 
three-dimensional Parker coronal heating model. HYPE- 
RION is a parallelized Fourier collocation-finite differ- 
ence code with third-order Runge-Kutta time discretiza- 
tion that solves the compressible MHD equations with 
DC field-aligned thermal conduction and radiation in- 
cluded. 

SETTING UP THE PROBLEM 
Governing equations 

We model the solar corona as a compressible, dissi- 
pative magnetofluid. The equations which govern such a 
system, written here in a dimensionless form, are: 
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where / = ^ [££yey + j(V x B) 2 — ^p 2 A(r)] . The 
system is closed by the equation of state, 

p = pT. (7) 

In the preceding equations the variables are defined in 
the following way: p(x,f) is the mass density, v(x,f) = 
(u,v,w) is the flow velocity, p(x,t) is the thermal pres- 
sure, A(x,f ) = (A x ,A y ,A z ) is the magnetic vector poten- 
tial, B(x,f) = (B x ,B y ,B z ) = V x A is the magnetic in- 
duction field expressed in terms of the associated Alfven 
velocity (B — > B/y/4npo), J = V x B is the electric cur- 
rent density, T(x,t) is the plasma temperature, £y = 
fl(djVj + d[Vj) — AV • \8jj is the viscous stress tensor, 
e ij = (dj v ' + diVj) is the strain tensor, and y is the adi- 
abatic ratio. The thermal conductivity (fc), magnetic re- 
sistivity (77), and shear viscosity (ju) are assumed to be 
constant and uniform, and Stokes relationship is assumed 
so the bulk viscosity A = (2/3) [1. The function g(z) de- 
fines the gravitational field strength: at t = we define 
p as poexp[— jgcos( j±)]. Assuming a uniform temper- 
ature we can determine the gravity as g = j3 i ^ . The 
function A(T) describes the temperature dependence of 
the radiation (A(T) = for T < 7b): 

A(T) = IzZ£ e («r-ro)/(TT ) _ t > Tq (8) 
Tq 

where Tq is the wall (photospheric) temperature, e = 2 
and T = 2 (fl]|2)). The important dimensionless numbers 
are: S v = PqVaLq/pl = viscous Lundquist number, S — 
VaLo/iI = Lundquist number, S r = VaLo/x = radiative 
Lundquist number (the new parameter % determines the 
strength of the radiation), j3 = Pq/Bq = pressure ratio 
at the wall, Pr = C p ji/K = Prandtl number, and A = 
Va/Vo = Alfven number. In these definitions, po is a 
characteristic density, Va is the vertical Alfven speed 
(used as the characteristic velocity to render velocities 
dimensionless), Lq is the vertical box length (= L z ), C p 
is the specific heat at constant pressure, C s is the free- 
stream sound speed, and Vb is the characteristic flow 
speed. Time (f) is measured in units of Alfven transit 
times (= Lq/Va)- 

Boundary Conditions and Forcing 

We solve the governing equations in a box of dimen- 
sions (L x ,L y ,L z ). The system has periodic boundary con- 
ditions in x and y, and line-tied boundary conditions in 
Z. To model a section of a coronal loop the system is 
threaded by a DC magnetic field in the z-direction (Bo). 

We then employ a simple, three-dimensional extension 
of the time-dependent forcing function used in the previ- 
ous studies iflOl 171. i.e., at the top and bottom walls we 
evolve a stream function: 

— J+/ 2 sin 2 (^— + -J (9) 



where fi(x,y) = Vol a' nm sin(k n x + k m y + Q m ). Val- 
ues for k are given by 3 < (k\ + kf n ) 1 < 4. At the top and 
bottom walls the magnetic vector potential is convected 
by the resulting flows. 

That is, the line-tied boundary conditions are: 

p(±L z /2) = p , 
pu(±L z /2) = -pody/dy, 
pv(±L z /2) = p dy/dx, 
pw(±L z /2) = 0, 
dA x /dt\ ±Lz/2 = vB Q , 
dA y /dt\ ±L _j2 = -uBo, 
B z (±L z /2) = B , 
T(±L z /2) = T . 
The enforcement of the boundary conditions is discussed 
in greater detail in 0. 

Numerics 

Equations 5 and 6 can be replaced by the magnetic 
vector potential equation: 

<9A 1 

— = vx Vx A+- Vx V x A (10) 
at S 

where A = V x B. Thus we solve numerically the equa- 
tions 1-3 and 7 together with equation 6. Space is dis- 
cretized in x and y with a Fourier collocation scheme 
with isotropic truncation dealiasing. Spatial deriva- 
tives are calculated in the appropriate transform space, 
and nonlinear product terms are advanced in configura- 
tion space. A second-order central difference technique 
H is used for the discretization in z. A staggered mesh 
also is employed in the z-direction[13|. In general, the 
fields that are defined at the z boundaries are advanced 
in time on the standard mesh. Other quantities of in- 
terest are defined and advanced in time on the stag- 
gered mesh. That is, on the standard mesh we look at 
p,pu, pv, pw, A x , A y , B z and T. Some derived fields 
such as (O x , (Oy, (O z , j x , and j y are also defined on 
the standard mesh. On the staggered mesh we look at 
A Xl B x , By, and j z . Note that for plotting purposes we in- 
terpolate these latter fields onto the standard mesh (at the 
boundaries an extrapolation is performed). 

A time-step splitting scheme is employed. All terms, 
with the exception of the vertical pressure gradient and 
the gravitation term, are discretized in time with a third- 
order Runge-Kutta scheme. The pressure step for the z- 
momentum is solved with a second-order Lax-Wendroff 
one-step central difference scheme. The vertical gravita- 
tion term is advanced using the forward Euler method. 

The code has been parallelized using MPI. A domain 
decomposition is employed in which the computational 
box is sliced up into x-y planes along the z direction. 
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FIGURE 1. Energies vs. time. Time is measured in units of 
axial Alfven crossing times L z /Va- 



FIGURE 2. Dissipation vs. time. Time is measured in units 
of axial Alfven crossing times L z /Va- 



NEW RESULTS 

In this section we report on the results of a prelim- 
inary numerical simulation of the model. This sim- 



ulation is run with L x = 2n,L y = 2n and L z 



87Z\ 



Other important parameters are g = 6.0,7 = 5/3,5 = 
S v = 80000,j3 = 0.001, B = 1.0, V = 0.01(^2/2), A = 
V A /V = 100\/2,po = 1.0,7b = 1.0,?* = 20.0,Pr = 
0.001, and S r = 0.0004. 

Temporal diagnostics 

To insure that we are obtaining good statistics, the sys- 
tem has to settle down into a steady state. Evidence for 
this is shown in Figure 1, which shows some of the im- 
portant energies as functions of time (time is expressed in 
units of Alfven transit times Lq/Va). As seen in our pre- 
vious RMHD simulations, the fluctuating magnetic and 

kinetic energies (e v = \ J^ z / 2 / 2 Sq } So" \v\ 2 dx dy dz and 

ej, = 2/^/2/0^ Jo* |b| 2 dx dy dz) settle down pretty 
quickly, with > e v (these quantities are also time in- 
termittent). Note, however, that the total internal energy 

Eint = :py J-/ 2 /2 Jo* Jo" P T dxdy dz takes much longer 
to level off in time. This reflects the fact that the system 
must heat up to attain the driven-dissipative steady state. 

The quantities shown in Figure 2 provide tem- 
poral information about the dissipation. Note that 
the radiation loss is a new quantity respect to 
previous simulations. Shown are the enstrophy 

1 rV 2 r L y r L * ni2 



J~ f-/ 2 /2 Jo 7 Jo* \ ®\ 2 dxdy dz, mean square electric 



current J = lo" ljl dx dy dz, and the total 

radiation losses D 



i! L f /2 !o y !o Lx P 2 MT)dxdydz 



as functions of time. The first two quantities behaves 
similarly to previous RMHD simulations, while the 
radiative losses settles on a similar timescale than the 
internal energy (Figure 1). 

Spatial diagnostics 

The following quantities are averaged over the perpen- 
dicular directions (x and y) and also over 1000 < t < 
2000. We look at these to determine the times averaged 
state of the system under unsteady heating. 

We first look at some of the quantities related 
to the dissipation of the system. Figure 3 shows 
some of the time averaged quantities as a func- 
tion of z, the direction of the large magnetic 
field Bo- Shown are the time averaged dissipa- 
tion intensity for the parallel vorticity Q z (z) =< 
[Jo Jo* 0) 2 (x 1 y,z)dxdy}^ > and also for the parallel 
electric current G z (z) = < [J^ Jq x j z 2 (x, y, z) dxdy] 2 > 
as well as the time averaged mean radiation rate 

A»(z) = 35: < Ip y Jo x p 2 A(T)(x,y,z)dxdy >. As in 
previous simulations, current density and vorticity 
are aligned to the dc magnetic field, so that their z- 
components are strongly dominant. Note that "< >" 
denotes the time averaging. The symmetry in z of 
these quantities indicates that we have averaged over a 
sufficient period of time. 

In Figure 4 we take a look at the time averaged ther- 
modynamic state for the unsteady heating case: shown 
are the time averaged mean mass density p m (z) =< 
Jo Jq x p (x, y, z)dxdy > and the time averaged mean tem- 
perature T m (z) =< J L> Jq x T(x,y,z)dxdy > as functions 
of z. 




The density profile is a result of the gravitational den- 
sity stratification. Figures 3 and 4 show feature typi- 
cal of the chromosphere-transition region-corona system, 
where density increases at lower heights, while tempera- 
ture increases in the high corona. Notice that most of the 
ohmic and viscous dissipation (Q z and G z ) takes place 
in the high corona, while radiation (D m ) origins mostly 
near the boundaries, where it is peaked (see Figure 3). 
This mostly results from the higher density values near 
the boundaries, as A(T) is multiplied by p 2 in the radia- 
tive term D m . 

DISCUSSION 

In this paper we have presented some preliminary results 
of our simulations of compressible DC coronal heating 
using our new HYPERION code. The inclusion of a 
thermal conductivity parallel to the DC magnetic field, 
coupled with a gravitational density stratification, gives 
rise to temperature and radiation features typical of a 
realistic coronal loop. 

This is an encouraging starting point to investigate the 
thermodynamical properties of a coronal loop threaded 
by a strong magnetic field whose footpoints are shuffled 
by photospheric motions. 
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